Rozpatrzmy równanie (zob. też moduł Nieliniowe rozwiązania okresowe w dwuwymiarowym układzie dynamicznym ):
\( \frac{d^2\,v}{d\,t^2}+\left(\beta v^2-2\alpha\right)\,\frac{d\,v}{d\,t}+v=0. \)
Równanie to można przedstawić w postaci układu dynamicznego
\( \left(\begin{array}{l} v \\w \end{array} \right)'= \left(\begin{array}{ll} 0 & 1 \\ -1 & 2\,\alpha \end{array} \right)\left(\begin{array}{l} v \\w \end{array} \right)+ \left(\begin{array}{l} 0 \\ -\beta\,v^2\,w \end{array} \right). \)
Rozwiązując równanie charakterystyczne
\( \det \left(\begin{array}{ll} -\lambda & 1 \\ -1 & 2\,\alpha -\lambda\end{array} \right)=\lambda^2-2\,\alpha\,\lambda+1=0, \)
otrzymamy:
\( \lambda_{\pm}(\alpha)=\alpha\pm\,i\,\sqrt{1-\alpha^2}=\alpha\,\pm\,i+O(|\alpha^2|). \)
Zatem, przy \( \alpha=0, \) macierz linearyzaji układu ( 2 ) w punkcie stacjonarnym \( (0,\,\,0) \) ma parę czysto urojonych wartości własnych \( \lambda_{\pm}(0)=\pm\,i \).
Przeanalizujemy teraz warunki powstania cyklu granicznego w otoczeniu tego punktu, kładąc \( \alpha =0 \). Obliczmy w tym celu macierz
której kolumny \( R=(R_1,\,R_2)^{tr} \), \( -I=-(I_1,\,I_2)^{tr} \) spełniają równanie
\( \left(\begin{array}{ll} 0 & 1 \\ -1 & 0 \end{array} \right) (R+i\,I)=iR-I. \)
Równanie to jest równoważne układowi
\( \left(\begin{array}{ll} 0 & 1 \\ -1 & 0 \end{array} \right) R=-I, \quad \left(\begin{array}{ll} 0 & 1 \\ -1 & 0\end{array} \right) I=R. \)
Zapiszmy pierwszy układ w postaci
Zatem uklad ( 3 )-( 6 ) jest spełniony gdy \( R_1=I_2 \), \( R_2=-I_1 \), przy czym \( I_1 \), \( I_2 \) mogą być wybrane w dowolny sposób (eliminujemy, oczywiście, przypadek \( I_1=I_2=0 \)).
Przyjmując, że \( I_1=-1 \), \( I_2=0 \) otrzymujemy następujący wzór:
\( P=(R,-I)=\left(\begin{array}{ll} 0 & 1 \\ 1 &0 \end{array} \right)=P^{-1}. \)
Mnożąc układ ( 1 ) z lewej strony przez \( P^{-1}=P \) oraz stosując zamianę zmiennych
otrzymujemy
(7)
\( \left(\begin{array}{l} x \\y \end{array} \right)'= \left(\begin{array}{ll} 0 & -1 \\ 1 & 0 \end{array} \right)\left(\begin{array}{l} x \\y \end{array} \right)+ \left(\begin{array}{l} -\beta xy^2 \\\,\,\,\,\, 0 \end{array} \right). \)
Zatem
(8)
\( f(x,\,y)=-\beta xy^2,\quad g(x,y)=0. \)
Wykorzystując wzór na obliczenie pierwszego indeksu Floqueta 1
\( \begin{array}{l} a=\frac{1}{16} \left[f_{xxx}+f_{xyy}+g_{xxy}+g_{yyy}\right]+ \\ +\frac{1}{16} \left[f_{xy}\left(f_{xx}+f_{yy} \right)-g_{xy}\left(g_{xx}+g_{yy} \right)-f_{xx}\,g_{xx}+f_{yy}\,g_{yy}\right], \end{array} \)
gdzie
\( f_{xx}=\frac{\partial^2\,f}{\partial\,x^2}(0,\,0), \quad f_{xy}=\frac{\partial^2\,f}{\partial\,x\,\partial\,y}(0,\,0) \quad i\,\,t.\,d. \)
otrzymujemy, że
\( a=\frac{1}{16}f_{xyy}=-\frac{\beta}{8}. \)
Rysunek 1: Portret fazowy układu uzyskany w eksperymentach numerycznych przy \(\alpha=0.2\), \(\beta=1\): pogrubiona linia czarna obrazuje rozwiązanie okresowe; lina czerwona oraz niebieska opisują trajektorie nawijające się na trajektorię okresową z obszaru wewnętrznego i zewnętrznego (odpowiednio).
Z powyższego wzoru wynika, że w przypadku, gdy \( \beta>0 \), układ posiada stabilne rozwiązanie okresowe przy małych \( \alpha>0 \). Przy \( \beta<0 \) oraz przy małych ujemnych wartościach \( \alpha, \) w układzie realizują się niestabilne rozwiązania okresowe.
Poniżej przedstawiono implementację w pakiecie "Mathematica" symulacji numerycznej równania Van der Pola. Na Rys. 1 przedstawiony jest portret fazowy ilustrujący przebieg rozwiązań w otoczeniu stabilnej trajektorii okresowej.
\( \begin{array}{l}Clear[\alpha ,\beta ,z] \\ z=1; \\ \beta =1;\\ \alpha =0.2;\\ \text{row1}=\partial _t V[t]==z*W[t] ;\\ \text{row2}=\partial _t W[t]==-z*( (\beta V[t]{}^{\wedge}2-2 \alpha )W[t]+V[t]);\\ \text{warV}=V[0]==0.01;\text{warW}=W[0]==0.0; \\ \text{rozw}=\text{NDSolve}[\{\text{row1},\text{row2},\text{warV},\text{warW}\},\{V,W\},\{t,0,100\},\\ \text{AccuracyGoal}\to 13, \text{PrecisionGoal}\to 13, \text{MaxSteps}\to 10{}^{\wedge}5];\\ \text{h30}=\text{ParametricPlot}[\{\text{Evaluate}[\{V[t],W[t]\}\text{/.}\text{rozw}]\},\{t,0,50\},\\ \text{PlotRange}\to \text{Full},\text{PlotStyle}\to \{\text{Red}\},\text{AxesLabel}\to \{V,W\}];\\ \text{warV1}=V[0]==1.75;\text{warW1}=W[0]==0.0;\\ \text{rozw2}=\text{NDSolve}[\{\text{row1},\text{row2},\text{warV1},\text{warW1}\},\{V,W\},\{t,0,100\},\\ \text{AccuracyGoal}\to 13, \text{PrecisionGoal}\to 13, \text{MaxSteps}\to 10{}^{\wedge}5];\\ \text{h31}=\text{ParametricPlot}[\{\text{Evaluate}[\{V[t],W[t]\}\text{/.}\text{rozw2}]\},\{t,0,30\},\\ \text{PlotRange}\to \text{Full},\text{PlotStyle}\to \{\text{Blue}\},\text{AxesLabel}\to \{V,W\}];\\ \text{warV2}=V[0]==1.27;\text{warW2}=W[0]==0.0;\\ \text{rozw2}=\text{NDSolve}[\{\text{row1},\text{row2},\text{warV2},\text{warW2}\},\{V,W\},\{t,0,100\},\\ \text{AccuracyGoal}\to 13, \text{PrecisionGoal}\to 13, \text{MaxSteps}\to 10{}^{\wedge}5];\\ \text{h32}=\text{ParametricPlot}[\{\text{Evaluate}[\{V[t],W[t]\}\text{/.}\text{rozw2}]\},\{t,0,20\},\\ \text{PlotRange}\to \text{Full},\text{PlotStyle}\to \{\text{Black},\text{Thick}\},\text{AxesLabel}\to \{V,W\}];\\ \text{fig2A}=\text{Show}[\text{h31},\text{h30},\text{h32}]\end{array} \)
Przypisy
1. Szczegóły wyprowadzenia tego wzoru można znaleźć w dopełnieniu do rozdziału 3.4 książki J. Guckenheimer and Ph. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Springer, NY, 1987
